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A method to calculate the bound states of three-atoms without resorting to 
an explicit partial wave decomposition is presented. The differential form of the 
Faddeev equations in the total angular momentum representation is used for this 
purpose. The method utilizes Cartesian coordinates combined with the tensor-trick 
preconditioning for large linear systems and Arnoldi's algorithm for eigenanalysis. 
As an example, we consider the He3 system in which the interatomic force has a 
very strong repulsive core that makes the three-body calculations with standard 
methods tedious and cumbersome requiring the inclusion of a large number of 
partial waves. The results obtained compare favorably with other results in the field. 
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I. INTRODUCTION 

In recent years the 4 He trimer has been the center of several theoretical investigations 
(see, for example, Refs. |2], |3[] and references therein). From all methods employed in 
these studies, the Faddeev equations method is perhaps the most attractive since it reduces 
the Schrodinger equation for three particles into a system of integral or differential equations 
which can be used to study bound states and scattering processes in a rigorous way. 

The differential form of the Faddeev equations has been proposed long ago by Noyes and 
Fiedeldey ||. Since then, these equations have been used in bound states calculations in 
Nuclear and Coulomb systems || |6], |7|, [| (| - just to mention a few references. The deriva- 
tion and discussion by Merkuriev of the asymptotic boundary condition of the differential 
Faddeev equations (DFE) ||10|| , paved the way to use the DFE not only in bound-state cal- 
culations but in investigations of scattering processes as well. For atomic systems, however, 
some peculiarities of the inter-atomic interactions resulted in a limited use of DFE in both 
bound and scattering calculations. 

One of the peculiarities is that the interaction among atoms often contains very strong 
repulsion at short distances which is difficult to handle numerically. In addition, it generates 
strong short range correlations implying that in atomic systems one should take into account 
many partial waves in order to achieve convergence. These problems can hardly be solved by 
using the computational power of modern computers and requires instead the development 
of appropriate numerical techniques. 

One such technique is the so-called tensor-trick proposed by the Groningen group [[IT], |T2| . 
Application of the method in Nuclear and Coulomb systems showed that high accuracy 
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calculations can be made using low computer power. Another method is that of Cartesian 



coordinates applied for the first time in three body Faddeev calculations in [13|. The latter 



method is suitable in describing the long-range behavior of weakly three-body bound states 
correctly. This feature is crucial in calculating excited states close to the two-body threshold. 
Yet another proposed method of solving the three-body problem is that of Ref. |T|| in which 
no explicit partial wave decomposition is required and only the total angular momentum 
of the system is used. Within this total angular momentum representation method the 
Faddeev operator has a very simple form that allows the construction of effective numerical 
schemes. 

In the present work we have combined the aforementioned methods and, moreover, we 
propose a new approach to overcome the difficulties arising from the short range repulsion. 
We applied the overall procedure in bound state calculations of the He 3 system using various 
potentials and compare our results with other results previously obtained by various groups. 

In Sec. II we describe the equations to be solved and formulate the appropriate boundary 
conditions. In Sec. Ill the numerical method is presented. Our results are given in Sec. 
IV while our conclusion are summarized in Sec. V. Some technical details for the field 
practitioners are shifted in the appendix. 



II. FORMALISM 

When considering a three-body system it is convenient to denote the two-body subsystems 
by a, (3 and 7 and introduce Jacobi coordinates which in the configuration space are defined 
by 



i^rj (r " _ri) (i) 

2m a (m l3 + m 7 ) \ / niprp + m 7 r 7 \ 
m a + mp + m>y J \ a nip + J 

The coordinates corresponding to other pairs can be obtained by cyclic permutations of the 
subscripts a, (3, and 7, their relation being 



X /3 \ _ ( *S/3a Spot 1 ' x < 



(2) 



where Sg a are coefficients which depend on the masses of the particles fllS] . For identical 
particles these coefficients are 

oil _ c22 _ 1 ol2 _ c21 

We assume that the Hamiltonian of the system involves only two-body interactions 

H = H + J2V l (x i ), (3) 

i 

where Hq = —A x — A y is the Hamiltonian of three free particles and Vi(xj) are the two-body 
potentials. 
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The wave function \1/ of the system can be expressed in terms of the three Faddeev 
components 

*(x,y) = y<) (4) 

i 

satisfying the Faddeev equations 

{-A x -A y + Vii^i) - E)^( Xi , y t ) = -Vfa) J2 ^-(xj, yj ) (5) 

where E is the energy of the system. In what follows we shall restrict ourselves, without 
loss of generality, to three identical bosons which is the case under consideration. In such a 
case the Faddeev components have the same functional form in their own coordinates and 
thus the system (|5|) is reduced to one equation only 

(H + V - E)<£> = -V(C + + C~)<S> . (6) 

In this equation Hq = —A x — A y , C + and C~ are the cyclic and anticyclic permutation 
operators respectively, $ is one of the Faddeev components written in its own coordinates. 

The potential energy of the system is invariant with respect to rotations. This makes it 
possible to separate out the degrees of freedom corresponding to rotations of the system. 
This can be achieved by expanding the Faddeev component $ in terms of eigenfunctions of 
the total angular momentum, i.e. Wigner functions D^ n (g) fl4|j , 

*(x,y)= E ^^^^ DU9). (7) 

L,m,n X V 

Here g G SO(3) are the coordinates describing collective angular motion of the system and 
(j) mn (x, y, z) are the projections of the Faddeev component in subspaces with fixed angular 
momentum. The components of the projections depend on the intrinsic coordinates 

s=l x l> l/=|y|> z = x,ye[0,oa), z 6 (-1,1) (8) 

xy 

which describe the intrinsic state of the cluster. Since we can fix the total angular momentum 
L of the system, the corresponding projection of the free Hamiltonian can be written as 

H£ = D L (g~ l )xyH —D L (g) (9) 

xy 

where D L (g) stands for a matrix constructed from the Wigner functions flTJ]]. In the case of 
zero total angular momentum the explicit expression for Hq reads 

d 2 d 2 1 1 d d 
H °° = ~^x~ 2 ~ df ~ ( ^ + " Z2) d~z ■ (10) 

The corresponding projection of Eq. (^) takes the form 

(H° + V(x) - E)4>°(x, y, z) = -V(x)P4>°(x, y, z) , (11) 



where 
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and x ± (x, y, z), y ± (x, y, z), and z ± (x, y, z) are 



+ / \ /"^2 3 9 V^3 \ 1 /9 

a {x,y,z) = (-x +-y ^f—xyz) 1 
y ± (x, y, z) = i^-x 2 + \y 2 ± ^-xyz) 1/2 



T ^ l y 2 xyz > ' (12) 
\/3 , a/3 2 1 



± — £ + — V xyz 

,±<„ v \ - 4 4 : 2 : 



^(x,?/,^) y±{x,y,z) 



Assuming that in each two-body subsystem only one bound state exists, we can write the 
asymptotic boundary conditions for the Faddeev component 0° as follows 

e _ fc 3 (x 2 +y2) l/2 

<f)°{x, y, z) ~ (p 2 (x) e~ k y y + A{x/y, z) ^ 2 + y2 y fi , (13) 

where <f2{x) denotes the wave function of the two-body bound state in the two-body 
subsystem, k y = y/E 2 — E 3 , k 3 = ^/—E 3l E 2 is the two-body bound state energy, and E 3 
the energy of the three-body system. The first term corresponds to virtual decay into a 
particle and a two-body bound system, usually denoted as 2+1, while the second term 
corresponds to a virtual decay with an amplitude A(x/y,z) into three single particles 
denoted as 1 + 1+1. The term corresponding to the latter configuration decreases much 
faster than in the 2+1. 

In the present work the term corresponding to 1+1 + 1 virtual decay is neglected and 
thus the asymptotic boundary conditions for the Faddeev component at sufficiently large 
distances R x and R y read 



d 

— In (p°(x,y,z) 



x=R x 



I — 9 
-k x =iJE 2 , —\n(f) (x,y,z) 

v dy 



= ~k y . (14) 

y=R y 



III. NUMERICAL PROCEDURE 

The first numerical calculation with the DFE were made by Laverne and Gignoux in 
the early seventies 0. Since then, many numerical methods were proposed. In the present 
work we shall employ: i) Quintique splines || together with an orthogonal collocation [ 16[ 



procedure which allows us to construct a linear system of equations corresponding to Eq. 
(|ll~l). In this way the bound state problem can be transformed to a generalized eigenvalue 



problem, ii) The tensor-trick UTTJ to solve the eigenvalue problem using the restarting Arnoldi 
algorithm [[17]. 

Any regular solution of Eq. ([Ll]) fulfilling the boundary conditions ( |14|) can be approxi- 
mated by an expansion in terms of the basic functions (see Appendix A) 

<fi°(x,y,z) =J2fiBi(x,y,z) (15) 

i 

where i stands for a multi-index {i x ,iy,i z } and Bi(x,y,z) = Bi x (x)B iy (y)Bi z (z). Following 
the procedure described in the Appendix A we obtain an equation for the coefficients ff 

(H F -EI)f = 0, (16) 
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where I is the unit matrix and Hp is the discrete analog of the Faddeev operator (see 
Appendix A) 

H F = H° + V(I + P). (17) 

When investigating nuclear systems with short-range potentials, such as the Malfliet-Tjon 
V (MT-V) potential the spectrum of the system (|I6|) can be calculated by direct 



application of the restarting Arnoldi or biorthogonal Lanczos algorithm. However, an 
additional regularization is usually needed to accelerate the convergence. For this, one 
may split the operator Hp in Eq. (|16|) into two parts, Hp = Hi + H 2 , where Hi, 
is such that the operator (Hi — EI) -1 can be explicitly constructed while H 2 can be 
considered as a remainder. The technique of explicit inversion of an operator Hi — EI 
is known as tensor-trick [TTfl and is briefly described in Appendix B for Cartesian coordinates. 



Following the standard procedure we write for the eigenvalue problem 

- (Hi -Ei)- x H 2 f = X(E)f . (18) 



where for physical solutions we have X(E) = 1. Eq. ( |18|) can be solved using the Arnoldi 
algorithm 



Although the above scheme can be applied effectively for a wide range of potentials, it is 
not satisfactory for interatomic potentials with a strong repulsive core. A typical example 
is the helium dimer potential which has a core with enormously strong repulsion and an 
extremely weak attractive well just enough to hold a two-body bound state. In what follows 
we give a brief analysis of the difficulties that arise and propose a recipe to overcome these 
difficulties. 

Since to reproduce the correct long-range behavior of the Faddeev components Cartesian 
coordinates are used, the natural choice of the operator Hi in (|1^) is 

Hi = H° + V . 

In this case the problem of calculating bound states is reduced to calculation of the largest 
eigenvalues of the operator 

L(z) = -(H + V - z)- l V(C + + C-) (19) 

To demonstrate how the presence of a short-range repulsion manifests itself in the spectral 
properties of the operators L(z), one must compare the spectra of potentials with and 
without a strong repulsive core. Two such potentials are the interatomic He-He TTYPT 
potential [|T^| and the nucleon-nucleon (MT-V) potential \TE\ the spectra of which are 



shown in Fig. [TJ. Two main features of the spectrum of the TTYPT potential should be 
pointed out, namely, that the spectrum contains numerous large negative eigenvalues and 
that for all values of the energy z, L(z) has a number of eigenvalues close to one. According 
to the estimations of the Arnoldi algorithm convergence rate (see Appendix C), these are 
the features that make the convergence unacceptably slow. 

In order to improve convergence we first look for the origin of the negative part of the 
spectrum. For this we consider the values such as Aj(zj) = —1 of the operator L(z). 
Evidently, they are simultaneously eigenvalues of the operator 

-^negative = Hq + V — VP . 
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Keeping in mind that the operator V stands for a potential with extremely strong short- 
range repulsion, we conclude that the term —VP produces extremely strong short-range 
attraction. Therefore, the operator -ff neg ative contains a large number of eigenvalues in the 
range [—00, E 2 } and their presence suppresses the convergence of the Arnoldi algorithm for 
values close to one. 

To eliminate the negative eigenvalues from the spectrum, we introduce, instead of L(z), 
the operator 

M{z) = {H + V + V m - z)- 1 (V m - VP) , (20) 



where V m corresponds to some strong short-range repulsive potential. It will be referred 
to as modifying potential. It can be shown that the eigenvalues \ii{z%i = 1 of the operator 
M(z) correspond to the eigenvalues z% of the original Faddeev operator and this does not 
depend upon the modifying potential. In contrast, the eigenvalues fii(z) = —1 correspond 
to the eigenvalues Zi of the equation 

~H + 2V m + V(l-P)-z t ]f = 0. (21) 

Obviously, these eigenvalues depend on the choice of the modifying potential V m a proper 
choice of which will eliminate all eigenvalues of the equation (|21| ) in the range z% G (—00, E%). 
Thus, according to Arnoldi algorithm convergence rate estimations, this feature of M{z) 
should dramatically improve the convergence to the maximal eigenvalue close to one. 

The convergence rate can be further improved if instead of the operator M a properly 
selected function of M is used in the calculations. One such function can be an even power of 
M(z) or the normalized Chebyshev polynomial of the second kind Uk(M(z))/Uk(l)- It allows 
an upward shifting of the lower bound of the spectrum and the increase of the separation 
between the physically interesting eigenvalues close to one and the rest part of the spectrum. 
According to the convergence rate estimations these features allow the reduction of the 
computation time and memory requirements considerably. 



IV. RESULTS 

We applied the aforementioned numerical procedure to study the spectrum of the He3 
trimer system. For this purpose, we employed the most recent He-He potentials, namely, 
the SAPT1, SAPT2, and HFD-B3-FC11 potentials |2T|, and compare the results with those 



previously obtained [fj| using the LM2M2 [2T| and TTYPT [B5| interactions, as well as with 
those obtained via other methods. 

In Table | we present the results for the trimer ground and excited states. We also 
present, in the same table, the results obtained for the dimer binding energy E2 calculated 
on the same grid used in three-body calculations as well as the exact results for the dimer 
E^ x . The difference between these values can be regarded as the lower bound for the error 
of our approach. 



In Tables |H] and [Hj we demonstrate the convergence, for the LM2M2 potential, of 



the calculated energies with respect to the number of grid points used. It is seen that 
convergence is achieved with a relatively low number of the z points (~ 12 points) while for 
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the other coordinates the number is much larger (~ 100 points). 

Our results for the LM2M2 potential are given in Table [TV] together with those obtained 
by Motovilov et. al. |22|, |23| via Faddeev-type equations (boundary condition model 
(BCM)) and by Nielsen et al. |IJ via the hyperspherical adiabatic approach. It is seen that 
an overall good agreement is achieved with these approaches. However, for the ground state 
a slightly better agreement is achieved with the results of Ref. |^3[ while for the excited 
state is with those of Ref. M. Comparing our results with the results of p| we should 
emphasize, that the angular basis used by Motovilov et. al. [E3 coincide with our spline 



basis for the simplest grid that can be used in the present method. Performing calculations 
with this simplified grid, we recovered all the digits of the ground state energy reported in 
p3| . This confirms the accuracy of their result and the suitability of the present method for 
bound state calculations. The better agreement for the excited state with the one obtained 
in Ref. |] can be attributed to the 2+1 contribution to the Faddeev component for the 
excited states. Taking into account this term in the BCM is very difficult, whereas the 
hyperspherical adiabatic approach of Nielsen et. al. |1] is more suitable for this purpose. 
However, the ground state energy of Jj] is about 1% less than our result. Consideration 
of the geometric properties of helium trimer can clarify the possible nature of the latter 
difference. 

The characteristic size of the bound states of the trimer can be estimated either by 
calculating (r) or (r 2 ) 1 / 2 . The results obtained for these radii are presented in Tables [VJ 
and [VI]. It is seen that they are approximately 10 times less than the radii of the dimer 
molecule. However, the size of the excited state has the same order of magnitude as that of 
the dimer. 

The wave function of the system can be easily obtained from the Faddeev component, 

,/ \ if x . (<P{x + ,y + ,z + ) (f)(x-,y-,z-)\ 

y, z) = <p{x, y, z)+xy — — + , 

\ x + y + x y J 

where x , y and z^ are defined by (fT2|). The most intuitive way to visualize the results of 
the calculations is to draw the one-particle density function defined as 



P (r 1 )=/dr 2 dr 3 |*(r 1 ,r 2 ,r3)P, 

V>(x(ri, r 2 , r 3 ), y(r h r 2 , r 3 ), z(n, r 2 , r 3 )) 



where 

#(ri,r 2 ,r 3 ) 

27Tx(r 1 ,r 2 ,r 3 )?/(r 1 ,r 2 ,r3) 

The functions x(ti, r 2 , r 3 ), y(rx, r 2 , r 3 ), and z(ri,r 2 ,r 3 ) are obtained from the Jacobi co- 
ordinates (HD and (§) and the wave function ip(x, y, z) is normalized to unity. Due to the 
symmetry properties of the system, the one-particle density function depends only on the 
coordinate r = |ri|. Taking into account the relation yi = v^ri we get 



p{r) = 47r2r2 J dxdz\ip(x,rV3,z)\' 



Omitting the integration over z, we obtain the conditional density function p(r, z) describing 
a spatial distribution for particle one when the other two particles are located along a fixed 
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axis. It is useful to plot this function in coordinates {r h r a ) such that 77 = rz is a projection 
of the position of particle 1 onto the axis connecting the other particles and 



r a = y^ r \^ ~ z r 

is a projection to the orthogonal axis. Three-dimensional plots of the function p((rf + 
r^) 1//2 , cosarctan ri/r a ) corresponding to the ground and excited states of the trimer calcu- 
lated with the LM2M2 potential are presented in Figs. and |3] respectively. The conditional 
density function of the ground state decreases in all directions in a similar way. The den- 
sity function of the excited state has two distinguishable maxima and exhibits the linear 
structure of the cluster. This structure has a simple physical explanation: The most prob- 
able positions of a particle in the excited state lie around the other two particles and when 
the latter particles are well separated the third one forms a dimer-like bound state with 
each of them. This interpretation agrees with the clusterisation coefficients presented in the 
Table |Vll|. These coefficients are calculated as a norm of the function f c defined by 



fc(y)= / dxdz<p{x,y,z)ip 2 {x) , (22) 



where (f2(x) is the dimer wave function. The values of \f c (y)\ , shown in the Table |VII| , 
demonstrate the dominating role of a two-body contribution to the trimer excited state 
whereas in the ground state this contribution is rather small. We could suppose that 
this dominating contribution of the cluster wave in the excited state ensured the fast 
convergence of the hyperspherical adiabatic expansion to the correct value, but in order 
to get the same order of accuracy for the ground state, possibly more basic functions should 
be taken into account. 

The advantage of using Faddeev equations over the Schrodinger one in bound-state cal- 



culations, is deduced from the results shown in Tables [VIII| and [La|. In the latter table we 
present the contribution of different angular states to the Faddeev component and to the 
wave function calculated as 

C n = \f n (x,y)\ 2 
fn{x,y) = J_ dzF(x,y,z)P n (z) , (23) 

where P n (z), n = 0,2,4, are the Legendre polynomials and F(x,y,z) is the Faddeev 
component or the wave function. The angular coefficients for the Faddeev component de- 
crease much faster than the wave function coefficients. This could explain the difference 
between the estimations of the trimer mean square radius reported in and the one pre- 
sented here. In the previous paper the same angular basis for the Faddeev component and 
the wave function has been used. However, it turned out that the contribution of the higher 
partial waves to the wave function is not negligible. This leads, for instance, to the change 
of the mean square radius of the excited state for the LM2M2 potential from 60.85 A to 
59.3 A. This observation agrees with the results of Nielsen et al. who reported the value 
60.86 A using essentially the same angular basis both for the wave function and for the 
Faddeev component. The Table |Vlll| also demonstrates that more angular functions should 



be taken into account in the ground state calculations. 
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V. CONCLUSIONS 

We present a method which can be used to perform accurate three-body bound states 
calculations. It is specially suited for systems where the interparticle forces have a 
strong repulsive core. The method is based on Faddeev equations in the total angular 
momentum representation and without any further partial wave decomposition. The 
equations are expressed in terms of Cartesian coordinates which are best suited to 
describe the long-range behavior of weakly three-body bound states. Combined with 
the tensor-trick preconditioning for large linear systems and the Arnoldi's algorithm for 
eigenanalysis it provided us accurate results for the the He3 trimer bound and excited states. 

Results obtained with the most recent realistic intermolecular forces (SAPT1, SAPT2, 
HFD-B3-FC11), as well as with earlier forces (LM2M2, TTYPT) indicate that only two 
bound states exist. The properties of these states are very different with the ground state 
being strongly bound while the binding energy for the excited state is comparable to that 
of the dimer. The latter implies that a dimmer cluster within the molecule is well formed. 

The sizes of these two states also differs much. The characteristic size of the ground state 
either estimated by (r) or (r 2 ) 1 / 2 is ~ 10 times less than the size of dimer molecule, but 
the size of the excited state has the same order of magnitude with that of the dimer. This 
estimation make it necessary to check for the absence of trimers in the experimental media 
during the measurement of dimer properties and vice versa. 

APPENDIX A: REDUCTION TO A LINEAR SYSTEM 

Suppose r t = {to, ti, . . . , t Nt } is a partition in the range [t ,t Nt ] in coordinate t and 
let 65,2(7*) be the space of quintique Hermite splines associated with this partition i.e. 
piecewise fifth order polynomial functions having two continuous derivatives in this range. 
The partition r t will be refereed to as the base grid. The set of basic functions in the space 
5*5,2 (t*) can be defined by the following conditions 

Sl{U) = l, d t Sl(U) = 0, #5?fe) = 

S?(ti)=0, flfc5?(ti) = l, #5?(*) = (Al) 

S*(U) = 0, d t S?(t t ) = 0, d?Sf(t t ) = l 

Vt ^ [imax(0,i-l),j ^min(i+l,AT t )] Sj(t) = 0. 

For simplicity we omit the subscript r^ t corresponding to any particular choice of a grid 
in each coordinate. Let the basic functions be arranged as follows: 

Bi(t) = Si(t), k = 0...N t , j = 1...3, i = 3k + j-l (A2) 

The solution of the Eq. (|ll|) must satisfy Eq. fll4]) and vanish on the planes x = and 
y — 0. To achieve this we modify the set of basic functions (|A1[) by constructing linear 
combinations Bi(t) which for t = x,y read 

S f (t)=Bi(t), z = 1...3JV t -l 

(A3) 

B 3Nt (t) = B 3Nt (t) - k t B 3Nt+1 {t) + k*B 3Nt+2 {t) 



10 



and for the coordinate z they coincide with the initial basic functions (|A2|) Bi(z) = Bi(z). 
We look for an approximate solution of Eq. (|TTD in the form of expansion in terms of the 
basic functions (|A3|) 

(f>°(x,y,z) = J2fiBi(x,y,z) (A4) 

i 

where % stands for the multi-index {i x , i y , i z } and Bi(x, y, z) = B ix (x)B iy (y)B iz (z). Since the 
basic functions (|A3|) fulfill the boundary conditions ([14]) these conditions are also satisfied 



by the approximate component 0°. 

To obtain a system of equations for the coefficients f\ the expansion (|A4|) should be 



substituted into Eq. (JlTj) with some subsequent projection procedure. A reasonable choice 
of such a procedure is the method of orthogonal collocations which do not require any 
integration. The highest order of approximation is achieved by a careful choice of collocation 



points for a given basic grid r pi. For the grids tn x , tm and one should construct 



the grids of collocation points r N c = {xl,x%,..., x c N c}, t N c = {y{, yf , . . . , y c N .}, and r NS = 
{zf, z%, . . . , z c N c} containing N£ = 3N X , N£ = 3A^, and N° = 3(N Z + 1) points respectively. 
Given these collocation grids the matrix elements of the operators involved in Eq. ([IT]) 
can be easily constructed. For example, the matrix elements of the identity operator are 
= Bi(x c jx ,y c jy ,z c jz ) with i = {i x , i y , i z } and j = {j x ,j y ,jz}- 
Substituting the expansion (Kl) into Eq. (O) followed by the orthogonal collocation 
procedure one gets a system of linear algebraic equations 

(F ° + V(I + P)- EI)f = . (A5) 

where Hq, V, I, and P are the matrices corresponding to the operators of Eq. (|TTJ) while f 
is the coefficient vector of the expansion ( |A4| ) . 

The dimension of this system is rather large and equals = N°NyN°. Therefore, we use 
this system of equations to find only its lowest generalized eigenvalues and eigenvectors. 



APPENDIX B: THE TENSOR-TRICK FOR CARTESIAN COORDINATES 

The inversion procedure of the operator (Hi — EI) is crucial for the performance of 
the algorithm and therefore a brief description of a tensor-trick algorithm for the DFE 
in Cartesian coordinates will be given. Furthermore, as the matrices involved are large, 
certain comments will be passed concerning the efficiency in computation. 

Consider the Hamiltonian H\ = + VI. The structure of this operator is 

T = T x (g)T y (g)T z 
% = D x ®I y ®I z + T x ®Dy®I z 

+ T X ®Iy®l Z (L x ®>ly + l x ® Ly) ®D Z + VI X ®Ty®T Z . 

The It are the matrices corresponding to the identity operators acting in coordinate t, It is 
the unit matrix, and D t contains the elements of the corresponding differential operators. 
The tensor structure of the operator H allows a diagonalisation of it without solving the 
spectral problem in the A^-dimensional space. To perform this diagonalisation one introduces 
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the following matrices B z , B x , B l y constructed from the generalized eigenvectors of the 
operators acting in different coordinates, 

B z : D Z B Z = L Z I Z B Z 

B x : (D x + VI X + kLJ x )B x = X%B X 

% : 05 y + kL y T y )B y = Y%B y 

where L z = diag{Zi, l 2 , . . . , l Nz }, X 1 = diag{E xl , E x2 , . . . , E xN J, and Y l = 
dia,g{E yl , E y2 , . . . , E yNy }, and E l tNt are the eigenvalues of the corresponding operators. One 

also defines the matrices B*, B* i} B* y constructed of the eigenvectors of the adjoined equa- 
tions. Expanding the expression B*(Hi — EI)B, with B and B* being given by 

N * „ 
B = U®(B x ®B y )®[B z ] t 
i=i 

N z 

b* = n©^®^)®^]*. 

and taking into account the relation B*B = \ x ® l y <g> l z , one can prove that the matrices B 
and B* are the desired matrices diagonalising the operator (Hi — EI). Once this operator 
is diagonalised it can then be easily inverted 

B*(Hi- EI)B = 6, 

(Hi - EI)- 1 = BG^B* 

G = diag{#m, g>n2, • • • , 9i x i y i z , • • • , Qn x n v n z } 

Ulxlylz ^Xlx ~ Xl y ^ 

Thus, in principle, the task of construction of the operator (Hi — EI)~ X is over. However, 
to use computer power effectively one must care about optimal implementation of all the 
operations involved in the computation of the product of the operator h~ l (E) = (Hi— EI)- 1 
and a vector u, the number of arithmetic operations should be minimized. To implement 
the operation effectively the structure of the operator h~ 1 (E) can be used 

h-^E) =BQ-\E)B*. 

The matrix G _1 (E) has explicit diagonal structure, and the calculation of the product G L u 
requires only 0(N x N y N z ) operations. Since however, the matrices B and B* are dense, 
and in principle the calculation of the products B'*'u requires 0((N x N y N z ) 2 ) arithmetic 
operations, that could essentially decrease the effectiveness of the method. However, taking 
into account the structure of the matrices B and B* one can improve the performance. 
Consider the procedure of tensor multiplication of two matrices by a vector u: 

v = B x <S> B y u , 

and use the identity 

B x ® B y — I <g) B y B x <g> I . 
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Obviously, even the calculation of the product B x ® B y u requires N x N y operations, the 
calculation of the products I® B y u and B x ® Tu require correspondingly only N x N y and 
N x Ny operations. Therefore a multiplication of a tensor product by a vector would cost 
N x N y (N x + N y ) operations. Similar relations can be written for a direct sum. Thus the 
correct implementation of a direct sum multiplication by a vector can reduce the computation 
costs N x N y N z /(N x + N y + N z ) times. 

The computation costs can be reduced even further. One of the possibilities is to abandon 
full diagonalisation of the operator h(z) and to modify the computation scheme. For this, 
let us construct the matrices diagonalising only the terms corresponding to the coordinates 
y and z 

N z 

B = n©(/«>5g® [B z ]i, 

1=1 

b* = n ©(/ ® ® • 

i=l 

In this case the matrix G(z) = B*h(z)B is band because of the fact that the operator is 
local and the basic functions have well localized support fl^3|). The multiplication L(z)u 
can be performed in two stages: 

i) Solve a system of equations with a band matrix G(E) 

G(z)u = B*(VP + V c S)u; 

ii) calculate v e I(z)u as 

v = Bu. 

Taking into account the tensor structure of the matrices B* and B this computation scheme 
requires only 0(N X N y N z (N y + N z )) operations. 



APPENDIX C: CONVERGENCE RATES FOR THE ARNOLDI ALGORITHM 

The key point of this algorithm is the repeated action of the operator under consideration 
on a sequence of vectors generating a Krylov subspace where the desired eigenvector lies. 
To use this algorithm to solve the eigenvalue problem ( |18|) one must be able to calculate the 
action of the operators (Hi — EI)^ 1 and H 2 on an arbitrary vector u. Thus the problems 
of inversion of the operator (Hi — EI) and storage of the matrix H 2 should be solved first. 
In Ref. [llj hyperspherical coordinates were used which lead to an effective storage scheme 



for the matrices H 2 = V(I + P). However, one can not include pair interaction into the 
invertible part (Hi — EI) using hyperspherical coordinates. On the other hand, if Cartesian 
coordinates are used then one has difficulties with the storage of the matrix P since these 
coordinates do not have the same symmetry properties as the operator P and this results to 
an irregular structure of P. In our case the situation is saved by the total angular momentum 
representation. In contrast to bispherical expansion used in PL the operator P is local and 



does not contain any terms requiring integration. As a result the matrix P is sparse and can 
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be either effectively stored using any standard storage scheme for sparse matrices |24| or not 



stored at all. In the last case the matrix elements are calculated "on fly" when required. 
The estimations of the convergence rate of the Arnoldi algorithm []17| are essential in 



understanding the problems that arise when the tensor-trick algorithm is applied. Suppose 
L is a linear operator in a finite-dimensional space with 

Lui = XiUi . 

Let /C m (v) be a Krylov subspace of dimension m constructed from a starting vector v of the 
operator L and V m be an orthogonal projector to JC m space. The convergence rate of the 
algorithm can be expressed in terms of the norm 

£ f) = ||(l-7> m H|| 

of the eigenvector u { projection onto the orthogonal completion to JC m . Suppose are the 
expansion coefficients of the starting vector v in terms of the eigenvectors Uj, 

n 

V = «iUj • 
i=l 

Then the following relation is fulfilled 

(TO) < £ (TO) 

fc i — SiS ) 

where & = — 1 + J2k=i \ a k\/\°£i\ an d e[ m ^ depends on the spectral properties of the operator 
L. For the eigenvalues lying in the external part of the spectrum some practically important 
estimations are known. Suppose that all the eigenvalues of L but Ai lye inside of an ellipse 
with the center c, largest semi-axis a, and focuses c ± e. It can then be shown that 
satisfy the following estimation ]17| 

, ^ Cm-l(-) 

,M |< _ e 



Cm-l( 



where C m (t) are Chebyshev polynomials of the first kind. 

Consider now the convergence rate of Arnoldi algorithm for the operator L(z) ([19|) for 
the value of z corresponding to the ground state of a three-body system. In this case the 
maximal eigenvalue Ai = 1. Having defined the separation distance between the largest 
eigenvalue and the rest part of the spectrum s = |1 — A2I, the distance between the maximal 
eigenvalues in the rest part d = |A 2 — A n | and neglecting the imaginary part of the spectrum, 
we obtain the estimation 

ie^i< 10^(1 +^)r i . 

Therefore the parameter managing the convergence rate is the ratio s/d of the separation 
to the size of the rest part of the spectrum. The convergence rate could be improved by 
enlarging the separation and by enhancing the remainder part of the spectrum localization. 
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FIGURES 




FIG. 1: Spectra of operators L(E) for potentials with (a) and without 
(b) a strong repulsive core 



'a 
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FIG. 2: Conditional one-particle density function of the He3 ground 
state, 77 and r a are given in A 




FIG. 3: Conditional one-particle density function of the He3 excited 
state, 77 and r a are given in A 
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TABLES 



TABLE I: Bound state energies for the helium dimer He2 and trimer He3: -Ef x and E 2 are the 
exact and calculated (using the same grid employed in three-body calculations) two-body bound 
states while E3 and E 3 are the ground and excited states for the trimer. 



Potential 



SAPT1|20| 
SAPT2j20| 
SAPT H 
LM2M2|21 



TTYPTp] 



Ef (mK) E 2 (mK) E 3 (K) E% (mK) 



-1.732405 
-1.815003 
-1.898390 
-1.303482 
-1.312262 



-1.7322 


-0.13382 


-2.790 


-1.8146 


-0.13516 


-2.888 


-1.8979 


-0.13637 


-2.986 


-1.304 


-0.12641 


-2.271 


-1.3121 


-0.12640 


-2.280 


-1.5873 


-0.13126 


-2.617 



TABLE II: Convergence for the He3 ground state energy with respect to the number of grid points 
for the LM2M2 potential. 

Grid E 3 (K) E 2 (mK) 



45x45x6 -0.12570 -1.3002 

60x60x6 -0.12583 -1.3032 

75x75x6 -0.12582 -1.3034 

90x90x6 -0.12582 -1.3035 

105x105x6 -0.12581 -1.3034 

105x105x9 -0.12636 -1.3034 

105x105x12 -0.12636 -1.3034 

105x105x15 -0.12640 -1.3034 

105x105x18 -0.12640 -1.3034 



TABLE III: Same as in Table || but for the excited state. 



Grid 


El (mK) E 2 (mK) 


45x45x6 


-2.2658 


-1.3011 


60x60x6 


-2.2649 


-1.3018 


75x75x6 


-2.2668 


-1.3031 


90x90x6 


-2.2670 


-1.3031 


105x105x6 


-2.2677 


-1.3037 


105x105x9 


-2.2712 


-1.3037 


105x105x12 


-2.2707 


-1.3037 


105x105x15 


-2.2707 


-1.3037 



TABLE IV: Comparison of the results obtained with the LM2M2 potential with other results 
the field. 



Observable This work 



E 3 , K -0.1264 -0.1252 -0.1259 

E%, mK -2.271 -2.269 -2.28 

< r 2 > 1 / 2 , A 6.32 6.24 

< r 2 > 1 / 2 , A 59.3 60.86 



TABLE V: The mean radius of He3 in A 

Potential Ground state of He3 Excited state of He3 He2 

SAPT1 5^47 47^5 45.60 

SAPT2 5.45 47.1 44.63 

SAPT 5.44 47.1 43.73 

HFD-B3-FC11 5.49 48.3 47.47 

LM2M2 5.55 49.9 52.00 

TTYPT 5.55 50.1 51.84 

TABLE VI: The mean square radius of He3 in A 

Potential Ground state Excited state He2 

SAPT1 6.22 56.4 61.89 

SAPT2 6.21 55.9 60.52 

SAPT 6.19 55.9 59.24 

HFD-B3-FC11 6.26 57.4 64.53 

LM2M2 6.32 59.3 70.93 

TTYPT 6.33 59.6 70.70 

TABLE VII: Contribution of cluster wave to the Faddeev component 

Potential ||/ c || 2 ||/ C *H 2 

SAPT1 0.2743 0.9442 

SAPT2 0.2787 0.9461 

SAPT 0.2832 0.9478 
HFD-B3-FC11 0.2661 0.9407 

LM2M2 0.2479 0.9319 

TTYPT 0.2487 0.9323 
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TABLE VIII: Contribution of different two-body angular states to the Faddeev component 

Ground state Excited state 

Potential S D G S D G 

SAPT1 0.9989914 0.0009974 0.0000108 0.9999951 0.0000049 0.0000001 

SAPT2 0.9989822 0.0010065 0.0000109 0.9999950 0.0000050 0.0000001 

SAPT 0.9989754 0.0010131 0.0000111 0.9999949 0.0000050 0.0000001 

HFD-B3-FC11 0.9990060 0.0009830 0.0000106 0.9999953 0.0000047 0.0000000 
LM2M2 0.9990393 0.0009500 0.0000103 0.9999956 0.0000043 0.0000000 

TTYPT 0.9990332 0.0009561 0.0000104 0.9999956 0.0000043 0.0000000 



TABLE IX: Contribution of different two-body angular states to the wave function 

Ground state Excited state 



Potential S D G SPG 

SAPT1 0.9521 0.0336 0.0089 0.855 0.099 0.031 

SAPT2 0.9520 0.0338 0.0090 0.854 0.100 0.031 

SAPT 0.9519 0.0339 0.0090 0.826 0.102 0.034 

HFD-B3-FC11 0.9525 0.0335 0.0089 0.837 0.100 0.033 

LM2M2 0.9530 0.0329 0.0089 0.873 0.093 0.025 

TTYPT 0.9527 0.0332 0.0094 0.862 0.094 0.029 



